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Abstract. In this work we study the heat transport in an XX Z spin-1/2 Heisenberg 
. chain with homogeneous magnetic field, incoherently driven out of equilibrium by 

reservoirs at the boundaries. We focus on the effect of bulk dephasing (energy- 
dissipative) processes in different parameter regimes of the system. The non- 
equilibrium steady state of the chain is obtained by simulating its evolution under 
the corresponding Lindblad master equation, using the time evolving block decimation 
method. In the absence of dephasing, the heat transport is ballistic for weak 
interactions, while being diffusive in the strongly-interacting regime, as evidenced by 
■ the heat-current scaling with the system size. When bulk dephasing takes place in the 

system, diffusive transport is induced in the weakly-interacting regime, with the heat 
current monotonically decreasing with the dephasing rate. In contrast, in the strongly- 



interacting regime, the heat current can be significantly enhanced by dephasing for 
, systems of small size. 



1. Introduction 



The transport properties of low- dimensional interacting quantum systems have been 
the object of intense research for several years. The development of powerful numerical 
methods have allowed important insights into the conduction regimes of these systems 
for different interaction strengths to be obtained [1, 2, 3, 4, 5]. In addition, the ground- 
breaking advances in the manipulation of cold atom systems have lead to the recent 
experimental observation of interesting transport phenomena, which strongly depend 
on the interaction between particles [6, 7, 8, 9, 10]. In spite of this vast theoretical 
and experimental progress, a complete picture of the transport properties of interacting 
quantum systems is still lacking. 

The analysis of transport phenomena becomes more involved when its coupling with 
the environment is considered. The study of the competition between this unavoidable 
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incoherent coupling and coherent dynamics, and its consequences on transport effects, 
has received much attention recently due to the observation of long-lived quantum 
coherences in light-harvesting complexes even at room temperature [11, 12, 13, 14, 15]. 
It has been observed that dephasing processes can enhance the transfer of particles 
[16, 17, 18, 19] and heat [20] in some non-interacting 2D networks, as well as in ID chains 
for non end-to-end conduction [21]; the optimal transport then occurs at an appropriate 
balance between coherent and incoherent effects. Recently we have shown that in 
a ID chain with end-to-end transport, a large enhancement of spin current emerges 
due to bulk dephasing, if strong interactions exist [22]$. Since dephasing processes 
take place in each site of the system, and their corresponding jump operators do not 
commute with the Hamiltonian, they not only randomize phase information between 
spins, but also dissipate energy. For interactions strong enough to produce separated 
bands of eigenstates of low mobility [23, 24, 25], dissipation induces transitions from such 
bands of large potential energy to scattering states with larger kinetic energy, thereby 
enhancing the spin current. This general mechanism of environment-assisted transport, 
described in Ref. [22], explains some recent experimental results, in which the presence 
of noise significantly enhances the expansion of a strongly-interacting bosonic gas in 
a disordered potential [7]. Other interesting dynamical consequences of the interplay 
between interactions and dissipation, e.g. interaction impeded decoherence [26, 27] and 
glass-like disorder-free states [28] have been lately reported. 

Experimentally, a commonly-used method to study transport properties of atomic 
quantum gases consists of releasing the atomic ensemble from its trap and observing 
its expansion under different conditions such as particle-particle interactions, disorder, 
noise, etc. [6, 7, 10]; several numerical simulations have been performed for similar 
situations, in both cold atomic gases [4] and spin systems [1, 2, 5]. A different method 
consists of analyzing the currents of systems driven out of equilibrium by asymmetric 
reservoirs at the boundaries [3, 29, 30, 31, 32, 33, 34, 35]. A very important step 
towards the experimental implementation of such configurations in atomic gases has 
been achieved recently, where the transport of ultracold fermions through a mesoscopic 
channel connecting two particle reservoirs has been realized [8, 9]. The incorporation 
of strong interactions by Feshbach resonances and of external noise by an appropriate 
modulation of an underlying optical lattice [7] is feasible given the recent advances in 
the control of ultracold atoms. 

An archetypical model to analyze the interplay of interactions and environmental effects 
on transport processes is given by the ID XXZ Hamiltonian, which is one of the 
simplest and most studied models of interacting quantum systems. This model can 
be experimentally realized in systems such as optical cavities [36], ultracold atoms in 
optical lattices [37] and some cuprate materials, with SrCu02 and Sr 2 Cu03 the most 
prominent realizations of the isotropic model [38, 39, 40, 41], and with (NO)Cu(NOs)3 
[42] and LiCuV0 4 [43] proposed to feature strong interactions. The commonly accepted 

| The model considered in Ref. [22] corresponds to a chain of interacting spinless fermions, which 
under a Jordan- Wigner transformation can be mapped to a spin chain with XXZ interactions. 
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picture of spin transport in this model, obtained from recent theoretical studies, 
corresponds to ballistic and diffusive behavior in the weakly- and strongly-interacting 
regimes, respectively [1, 2, 3, 29, 44]. Theoretical analysis, based on the integrability 
of the XX Z model [23] and the conservation of the heat current [45], has led to 
the conclusion that for any interaction strength between spin excitations the thermal 
Drude weight of the system is nonzero. Consequently the heat conductivity is infinite, 
i.e., the heat transport is ballistic [45]. This result agrees with several analytical 
[46, 47, 48, 49, 50, 51, 52, 53, 54] and numerical [5, 55] studies, and has been supported 
experimentally by the measurement of large heat conductivities in cuprates [56, 57, 58]. 
Interestingly, the recent theoretical studies on non-equilibrium XXZ chains using 
boundary reservoirs suggest that spin and heat transport properties differ from this 
generally accepted picture. A transition from ballistic to diffusive heat transport for 
strong interactions between spin excitations was observed in a chain with a strong 
magnetic field [30]. In a similar setup with lower fields, ballistic heat and spin transport 
were observed for the non-interacting regime only, arguing that finite interactions induce 
diffusive transport [31]. Therefore more effort is still required to fully understand the 
transport properties in these driven systems. 

To help establish a complete picture of transport phenomena in quantum systems where 
interactions and dephasing have a prominent role, in the present work we study the 
heat current through an XXZ chain driven out of equilibrium by unequal reservoirs 
at the boundaries, and with dephasing processes in the bulk. Initially we present the 
model to be analyzed and the definition of the local heat current in Section 2. We 
then study the heat transport properties of the XXZ spin chain for several parameter 
regimes. In Section 3 we characterize the nature of heat transport through the system 
in a homogeneous magnetic field without dephasing. We find that heat propagates 
ballistically in the weakly-interacting regime, while for strong interactions, negative 
differential conductivity (NDC) emerges, with diffusive behavior at weak driving only 
and an insulating state at large driving. Then we observe the effects of bulk dephasing 
on heat transport in the different interaction regimes. In Section 4 we consider the 
case of weak interactions. By analogy with Refs. [59, 60] we show that dephasing 
induces a non-equilibrium phase transition from ballistic to diffusive regimes, with the 
heat transport monotonically degraded by dephasing. For strong interactions, studied 
in Section 4, dephasing degrades the NDC effect, leading to a diffusive behavior of the 
system for all driving strengths. In addition, the heat current can be enhanced by 
moderate dephasing rates for small systems. 
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2. Model of heat transport in a ID Heisenberg chain 

2.1. N on- equilibrium XX Z spin chain 

We consider the nearest-neighbor interacting spin-| XX Z spin chain, in the presence 
of a homogeneous magnetic field. The Hamiltonian is given by (taking h — 1) 

JV-1 N 

H = ^ 5>K« + + ^U) + * £ (!) 

i=i j=i 

with erf (a = x,y, z) the Pauli matrices at lattice site i, N the number of sites, r the 
exchange coupling between nearest neighbors (we set the energy scale by taking r = 1 
from here on), B the magnetic field and A the anisotropy parameter. The (t|(j| +1 terms 
act as an interaction penalizing the alignment of neighboring spins, while the in-plane 
terms cause spin-flip "hopping" processes. The ground-state phase diagram of this 
model is well known: in the absence of magnetic field, A < — 1 and A > 1 correspond 
to ferro- and antiferromagnetic (gapped) phases, respectively, and — 1 < A < 1 to a 
gapless state. A finite B shifts the critical points to a larger interaction strength |A| if 
A > 0, and to a weaker strength |A| if A < 0; see Ref. [61] for a detailed description. 
In the present work we will only consider interactions A > 0, and magnetic field B > 0. 
To model a non-equilibrium configuration we introduce Markovian reservoirs. The 
dynamics of the system is captured by a Lindblad quantum master equation [62] 

^ = -i[H,p] + £(p), (2) 

where p is the density matrix of the chain and C(p) is the dissipator describing its 
coupling to the environment. In Lindblad form, 

c( P ) = = E ( L *p L t - \i L i L ^ p}) > ( 3 ) 

k k ^ ' 

where Lj, are the jump operators corresponding to the coupling of the system to the 
environment, and {.,.} is the anti-commutator. We consider three different couplings 
C(p) = + Ai(p) + Ai(p), described below. 

The dissipators £l(p) an d £r(p) represent the boundary driving of the system to a 
non-equilibrium configuration, and correspond to the jump operators [29] 

L£ tR = y/T(l±f)/2*t,N L^ R =^r(l T f)/2a- N , (4) 

with = \{p X j ± i a j)- These operators create and annihilate spin excitations at the 
boundary sites of the system, with / the driving strength (0 < / < 1) and V the coupling 
strength to the spin reservoirs (we take T = 1). The driving strength / controls the non- 
equilibrium forcing, and in isolation would impose unequal average spin polarizations 
in the boundaries, (af) = f and {o~ z N ) — —f [60]. A zero bias / = leads to a steady 
state p = t/2 N [63], indicating that the system is formally at infinite temperature; 
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Figure 1. Scheme of a spin chain driven out of equilibrium by the action of reservoirs 
at the edges, with coupling strength T and driving parameter /. The dephasing 
processes have a rate 7. The light arrows indicate the homogeneous flow of spin 
excitations. The dark arrows represent the heat current, which varies through the 
chain due to the energy dissipation; the amplitudes of the latter are based on the 
results shown in Fig. 6(b). 



no spin nor heat current exists in this case since spin-flip excitations are created and 
annihilated at the same rate on both boundaries. At finite driving / a net spin current 
flow is established due to the asymmetry between the two edges. At maximum bias 
/ = 1 the boundary jump operators reduce to and a N . 

The bulk dephasing processes are described by the dissipator £d(p) = X^Li^d 1 G°)> 
with a homogeneous dephasing rate 7 and jump operator at site j 

L A 3 = s/iv*. (5) 

A schematic representation of the system, including the different incoherent processes, 
is shown in Fig. 1. 

The interplay between coherent and incoherent phenomena leads to a non-equilibrium 
steady state (NESS) p ss , which satisfies dp ss /dt = 0. This state is thus identified by the 
stationary values of its spin and heat currents. Since the Hilbert space containing 
p ss grows exponentially with the size of the system (as 4^), it is not possible to 
obtain the exact steady state for large N numerically. In addition, a simple analytical 
approximation to the steady state of a large system is only feasible in the non-interacting 
regime [44, 60]. To obtain p ss , we then directly simulate Eq. (2) and take the long 
time limit applying the time evolving block decimation (TEBD) method [64, 65] to a 
matrix product operator description of p(t). Our implementation of the algorithm is 
based on the open source Tensor Network Theory (TNT) library [66], and allows us 
to analyze efficiently the transport response of strongly-interacting systems larger than 
those considered in previous studies [30, 31]. 

2.2. Alternative driving schemes 

In this work we study the heat transport induced by a magnetization imbalance at the 
boundaries of a spin chain, imposed by the Lindblad driving processes of Eq. (4). This 
means that we actually observe the magnetothermal behavior of the system, instead 
of the direct response to an energy imbalance enforced at the boundaries (see also 
[67]). To analyze the latter situation, different driving mechanisms have been proposed. 
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One possibility consists on using two-spin driving operators [3], designed to induce a 
Gibbs state of determined temperature when applied to an isolated pair of interacting 
spins (analogous to the single-spin operators of Eq. (4) which induce a state with 
magnetization / on an isolated spin). Another recently proposed configuration consists 
of connecting the chain to finite meso-reservoirs at the boundaries, driven to a state 
of fixed temperature and chemical potential by a Lindblad-type coupling to super- 
reservoirs [55]. 

These alternative and somewhat complicated driving processes would allow a complete 
study of the heat current directly generated by an energy imbalance imposed at the 
boundaries of the system, as well as of its magnetothermal response, for different 
interacting regimes. Here we instead consider the much simpler and intuitive driving 
scheme of Eq. (4), which still allows the observation of interesting phenomena in various 
parameter regimes. The transport properties of strongly-correlated systems under the 
alternative driving schemes are the subject of current research [68]. 

2.3. Spin and heat currents 

In our previous work [22], we analyzed the spin transport properties of the system 
described above, by computing the spin current and local magnetization. The spin 
current operator follows from the continuity equation, and is given by 

^ = 2(« 1 -af< +1 ), i = l,...,JV-l, (6) 

where its expectation value is spatially homogeneous in the steady state, so (Jf) ss = 
Tr (p ss jf) = J s for all i§. The spin current J s and the magnetization profile indicate 
different spin transport regimes. The ballistic regime, corresponding to A < 1 and 7 = 0, 
is characterized by a nearly-flat magnetization profile, and by a current proportional to 
the driving and independent of the size of the system, i.e. J s oc / [60, 59]. In the 
diffusive regime, corresponding to A > 1, weak driving and 7 = 0, as well as to any A 
and 7 > [22, 59], the system satisfies the diffusion equation 

J S = k s VS, (7) 

with k s the spin conductivity and VS = ({erf) — (a z N ))/(N — 1) the magnetization 
gradient. This transport regime is characterized by a linear magnetization profile and 
a spin current scaling J s oc f/N [22, 29]. 

In the present paper we focus on heat transport. Due to the energy dissipation caused 
by dephasing, the heat current is not homogeneous in the steady state, but different at 
every site. We denote by Jf the heat current operator at site i; to obtain its expectation 

§ In the present work we are interested in the transport properties in the steady state of the system. 
Hereafter all the expectation values are calculated for such a state, and for simplicity we will drop the 
subindex ss indicating steady state expectation values. Nevertheless the definitions of heat currents 
presented in Section 2.3 remain valid in the transient regime. 



Heat transport in the XX Z spin chain 

value, we first write the Hamiltonian in the form 



N-l N-l 



H = ^£,,,+1 = ^2(h iji+ i + b iA+1 ), (8) 

i=l i=l 

where h iji+ i and b iji+1 correspond to the XXZ coupling and magnetic field terms, 
respectively, and are given by 

h hl+l = a*a* +1 + a y t a v l+1 + Aa?<7? +1 , & M+ i = -[<r?(l + 8 iA ) + <r? +1 (l + <S i+ i,iv)], (9) 
where 5 denotes the Kronecker delta. 

To obtain the heat current, we calculate the rate of change of the local energy density 
£i,i+i. From the master equation (2), we have 

d(e iti+1 ) _ / <9p 



9t 



T^e M+1 |^ = -iTr(e^ +1 [H,p]) + Tr(£ M+1 £ L (p)) 



AT (10) 

+ Tr(^ +1 £ R p)) + ^Tr(^ +1 4 j) (p)), 

j'=i 

where i = 1, . . . , N — 1. In the steady state, the previous equation is reduced to 

^±i> = i([#,£ M+1 ]> + Tr(£ 1 , 2 £ L (p))^, 1 + Tr(£ JV _ 1;Ar £ R (p))^ +liiV 

+ Tr^+^p)) + Tr(e,, +1 4 +1) (p)) = 0. 

Next we consider the continuity equation at each pair of neighboring sites in the steady 
state 

d{e iii+ i) H 9(£ i)i+ i) 



-V(Jf ) + 



= 0, (12) 

Env 



9t x 1 ' <9t 

where V(J 4 H ) = (J^x) — (Jf) is the gradient of the heat current, and (d{e i:i+1 ) / dt)\ En v 
is the rate of change of the local energy density due to the coupling of the system with 
the environment. By comparing Eqs. (11) and (12), we identify 

i{[H,e i>i+1 \) = -V(Jf) = (Jf)-{Jf +l ) (J?) = *([£,_!„ £ M+1 ]>, 2 < i < N-l. 

(13) 

Note that this definition of the heat current does not apply to the boundaries of the 
chain. For such cases, the heat current results from the action of the left and right 
reservoirs. If we take 

(J?) = Tr(e li2 £ L (p)) - (jj=> = Tr(£^_ liJV £ R (p)), (14) 

then we have a consistent definition of the heat current over the entire system, since the 
continuity equation for any pair of nearest neighbors gives 

<■#!> - = Tr(£ M+1 £«(p)) + Tr(£ M+1 £f 1} (p)) = -4 7 <« + i + afaf +1 ). (15) 
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As expected, the difference between heat currents at neighboring sites arises solely from 
the energy dissipation due to dephasing processes at those sites. 

In order to understand the features of the heat transport we separate the heat current 
in Eq. (13) into two contributions: that arising from XXZ interactions only, given at 
site % (2 < % < N - 1) by 

(</f XZ > = i([hi-i,i, K l+1 }) = 2((al 1 aX +1 - , . 

and the heat current carried by the net flow of spin excitations, corresponding to 

( J ? ) = *([e»-M, ~ i{[hi-i,i, Ki + i}) 

= B{(ot_ ia ? - at,*?) + (*?a! +1 - afaf +1 )> = |((J^ + J?)> (17) 

We define similar contributions at the boundaries. For site i — 1, 

(J xxz > = Tr(/i li2 £ L (p)), (Jf) = Tr(6 li2 £ L (p)), (18) 

where we obtain 

(j xxz > = -r«M + a( ( t 1 V|))/2 + r/A(o5>, (Jf) = rB(/ - K». (19) 

The corresponding definitions for site are 

(J xxz ) =r((ft JV _ 1 ^) + AK. 1 ^))/2 + r/AK_ 1 > ) (j»> = rs(/ + K». (20) 

The total heat current operator at site i is then J t H = J xxz + J?. Since the spin 
current of the steady state of the system is homogeneous, so is Jf\ we thus define 
J B = (Jf) = BJ S . The spatial variation of the heat current arising from dephasing 
processes is then entirely contained within (J 4 XXZ ). 

Note from Eqs. (19) and (20) that the spin current J s and the homogeneous component 
of the heat current J B can be directly obtained from the magnetization at the 
boundaries. This was previously found for the case A = [60]. Now we observe 
that the same applies for all interaction strengths A, and that J s and J B are given by 
the deviation of the magnetization of the boundary spins from the one that would be 
imposed by the left and right reservoirs to isolated spins (±/). 

We now study the total heat current and its components in the non-equilibrium setup 
described in Section 2.1, and observe its behavior in the weakly (A < 1) and strongly 
(A > 1) interacting regimes, in the absence and presence of dephasing processes. 



3. Heat transport in driven XXZ spin chains without dephasing 

Initially we consider a system without dephasing processes. We note that if the magnetic 
field is absent (B = 0), there is no net heat flow through the chain, in spite of the 
existence of a finite spin current. This can be understood qualitatively by considering 
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Figure 2. Schematic representation of heat transport in the absence of dephasing. 
(a) If B = 0, the heat current of spins up injected at site 1 (dark arrows) cancels with 
that of spins down injected at site N (light arrows), and the net current is zero, (b) If 
B > 0, the first current becomes larger than the second, and a net current emerges. 



the strongly-driven case / = 1, in which only spins up are injected at site 1 and only 
spins down are injected at site iV (the argument is easily extended to / < 1). At 
each site of the chain, the total heat current is the superposition of two contributions: 
one current flowing from right to left, carrying energy associated with the spins down 
injected at site N, and other flowing from left to right, carrying energy associated with 
the spins up injected at site 1. Without any energetic asymmetry in the system, the 
contributions cancel with each other since there is no preference for any z direction, as 
illustrated in Fig. 2(a). If a homogenous magnetic field is turned on, an asymmetry is 
established, and a total heat current emerges. If B > (as in the present work), the 
left-to-right flow carries energy of spins aligned with the field, so it becomes larger than 
the right-to- left flow. The net heat current is thus positive, as schematized in Fig. 2(b). 
The vanishing of the heat current at 7 = and B = can be formally demonstrated by 
considering the symmetries of the Lindblad master equation (2) with XX Z Hamiltonian 
and the boundary driving of Eq. (4), as discussed in Ref. [67]. In this case, the 
NESS of the system is invariant under the transformation U = Q a R, where Q a = 
o~x <E> o~% Cg> • • • <g> a%, a = x,y, and R is a reflection operator, so R(Ai Cg> B 2 <E> • • • <E> Cn) = 
(Ci ® • • • <g> fijv-i <g) A N )R [67]. This means that 



P 



UpU^ = Q a RpRQ a . (21) 



Since J 4 XXZ changes sign when transformed by U Jf^U = — J 4 XXZ ), the symmetry 
of Eq. (21) leads to the following property [67] 

(J xxz } = Tr(J xxz p) = ^(jxxz^t) = ^uUf^Up) 

= -Tr(J xxz p) = -(J xxz ) (J xxz ) = 0. 1 j 

So in the absence of dephasing and magnetic field, both components of the heat current 
vanish, and (J? ) = for all i. 

We now discuss the results of heat transport for B > 0. The first important observation 
is that the spin current remains unmodified if a homogeneous magnetic field is present 
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in the system. This occurs because the effect of the magnetic field is to vary the relative 
energies of the sectors of fixed spin quantum number in the spectrum of the model||, 
while leaving the internal structure of each sector unchanged. This does not affect the 
spin transport through the chain, which is governed by processes that conserve the 
number of spin excitations. The second important point is that similarly to the spin 
current, (J xxz ) is independent of the magnetic field. This can be seen by rewriting Eq. 
(16) in terms of a + and a~\ for example, the first term is 

°i-i<rZ°i+i - °i-i<rZ°i+i = 2*(<7~ i«i " °ti° Z i°7+i)- (23) 

The resulting terms are spin-conserving, so they only involve processes within each spin 
sector separately. A magnetic field does not affect expectation values of spin-conserving 
operators like J xxz . Therefore for any homogeneous magnetic field in the absence of 
dephasing, (J xxz ) = 0, and the heat current becomes homogeneous through the chain 
and possesses all the features of the spin current, since (Jf 1 ) = J R = J B = BJ S . Note 
that this is reflected in the absence of spatial imbalance of the energy density inducing 
J xxz , i.e., (hij+i) is homogeneous through the system. 

It is important to emphasize that the symmetry of Eq. (21) specifically applies to the 
Lindblad driving scheme of Eq. (4). Different simple driving configurations might lead to 
similar symmetries, resulting in vanishing particle or energy currents [67]. Nevertheless, 
for other types of driving, these symmetries do not apply, leading to contributions of 
both J B and (J xxz ) to the total current even in the dephasing- free case, and thus to 
different transport regimes than those observed in this work. 

We now study the transport properties of the dephasing-free chain at weak driving and 
different interaction strengths A. Observing the scaling of the current with the size of 
the system, we can identify different transport regimes. As shown in Fig. 3(a) for A < 1, 
the heat current remains constant as N increases, indicating ballistic transport. This 
picture is reinforced by the corresponding nearly flat energy profiles presented in Fig. 
3(c). A different situation is observed for A > 1; as seen in Fig. 3(b), the heat currents 
decrease with the size of the system, and as shown in Fig. 3(c), the corresponding 
energy profiles have a constant finite gradient away from the boundaries of the chain. 
This behavior is characteristic of diffusive transport, in which the current is given by a 
diffusion equation 

AE 

J H = Kll VE with VE = — . (24) 

N -4 v ' 

This is the discrete version of Fourier's law, with k h the heat conductivity, VE = 

(£i,i+i) — = (bi,i+i) ~ (h-i,i) the constant energy gradient in the bulk and 

AE = (en-2,n-i) — (£2,3) = (on-2,n~i) — (^2,3) the energy difference across the chain 

(excluding the boundaries). This equation indicates that the heat current decays with 

the size of the system as J H /AE ~ TV -1 , which is found to agree very well with the 

|| Each quantum number is an eigenvalue of the total spin in z direction ^ of, which is conserved in 
the XX Z model. 
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Figure 3. Heat transport properties of a chain with homogeneous magnetic field and 
no dephasing as a function of N, with driving parameters / = 0.1 and T = 1, for 
different interaction strengths, (a) Heat current for A = 0.9 (B = 0.5) and A = 0.5 
(B = 1.0); the lines are guides to the eye. (b) Log- log plot of the heat current for 
A > 1. The lines correspond to the fit J n /AE = k/(N - A) a . For A = 2.0, a = 0.98 
and k = 0.85; for A = 1.2, a = 1.02 and k — 2.95. (c) Corresponding energy profiles 
for N — 40, where A increases in the direction of the arrow. 

current scaling shown in Fig. 3(b). Thus for the interactions A > 1 considered, the 
diffusion equation (24) is satisfied. In addition, we observe that 

J H = B J S = « H ((6 i)i+1 ) - (&<_!,<)) = y« H («fi> - (trti)) = B K}1 VS. (25) 

This means that the heat and spin conductivities are equal, i.e. kh = K s- 
The results shown in Fig. 3 indicate that the critical point separating ballistic and 
diffusive transport regimes is very close to the isotropic point. In fact, since the heat 
current of the present configuration is entirely given by the magnetic field and the 
spin current, the critical point corresponds to A = 1, which separates the ballistic 
and diffusive regimes of spin transport [1, 2, 29]. This result contrasts with previous 
analysis of heat transport through chains driven by external reservoirs, which indicate 
ballistic transport for the non-interacting case only [31] and for interactions A < 1.6 
[30]. It also differs from the commonly-accepted result of finite thermal Drude weight 
and ballistic heat transport for all interactions strengths A [5, 45], supported by 
integrability arguments. Nevertheless, the introduction of external reservoirs can break 
the integrability of the system in the most general case^f [31, 32], so it is not expected that 
the picture of ballistic transport for all A maintains its validity in the non-equilibrium 
configuration considered in the present work. 

It is also important to note that the ground state phase diagram of the system does 
not determine the nature of the (spin or heat) transport response seen here, since 

% Up to now, only the non-interacting case A = of the boundary-driven configuration considered 
here has been shown to be integrable [69, 70]; it is not known whether interacting cases are integrable 
too. 
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Figure 4. Transport properties of the strongly- interacting case A = 1.5 and driving 
parameters / = 0.1, 0.2, . . . , 1.0, in the absence of dephasing; B = 0.5, T — 1 and 
N = 16. (a). Heat current, (b) Profiles of energy, (c) Magnetization. In (b) and (c), 
the increase of / is indicated by the arrows. 

the equilibrium critical points strongly depend on the magnetic field [61], while the 
interaction separating the different transport regimes is field-independent (A = 1). 
Thus, the transport properties of the chain are determined by the structure of its 
eigenspectrum, which qualitatively changes from being continuous for A < 1 to 
consisting of separated bands for A > 1 [23]. 

Now we focus on the strongly-interacting regime, and obtain the heat current of the 
system for several drivings. Since the convergence to the NESS is slow at large driving 
/ ~ 1 [22, 29], we consider systems of small size, namely N = 16. Nevertheless, test 
calculations for larger systems present a similar behavior. The results are shown in 
Fig. 4. For weak driving, the heat current increases linearly with / (see Fig. 4(a)), 
decreases with the size of the system (as shown in Fig. 3(a) for / = 0.1), and the 
energy profiles have a homogeneous gradient in the bulk (see Fig. 4(b)), indicating 
diffusive transport. As the driving increases beyond / ~ 0.5, the current starts to 
decrease. This large imbalance between injection and ejection of spin excitations at the 
boundaries favors the occupation of flat bands of bound states with low mobility (which 
only exist for A > 1), making it harder for the excitations to move through the chain 
and thus diminishing the current; this mechanism is described in Ref. [22]. At the 
strongest driving / = 1 the heat current vanishes, corresponding to insulating behavior 
(in fact it is exponentially suppressed with N, a tendency proved for the spin current 
[71]). This phenomenon of decreasing transport with increasing driving, non-existent for 
weak interactions, is known as negative differential conductivity (NDC), and has been 
previously observed in the heat transport of XY chains with B > [70] and in the spin 
current through an XX Z chain with B = [22, 29]. In the latter case, the insulating 
behavior at very strong driving results from the appearance of ferromagnetic domains 
of almost entirely-polarized spins at the edges of the chain, which strongly inhibit spin 
flips. Since the presence of a homogeneous magnetic field does not modify the spin 
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Figure 5. Effect of dephasing on the counter-propagating heat currents coming from 
each boundary. When reaching each site, the currents have been attenuated by a 
different amount of dephasing processes, so their local magnitudes are distinct. 



transport features, the existence of ferromagnetic domains is maintained, as shown in 
Fig. 4(c), leading to the absence of both spin and heat currents at / = 1. 

4. Effects of dephasing on heat transport for weak interactions 

In this section we discuss how dephasing processes affect the heat transport properties in 
the weakly-interacting regime. It is illustrative to study first the system in the absence 
of magnetic field (B = 0), where the local heat current is entirely given by (</ xxz ) 
(since J B = 0). We consider again the picture of two currents flowing in opposite 
directions, associated with spin excitations injected at each boundary. In this case, the 
currents do not cancel with each other. Instead, they reach each site of the system 
with unequal magnitudes, depending on the attenuation they have experienced due to 
dephasing processes along the different distances they have covered, as depicted in Fig. 

5. As a result, the net current in the left half on the chain is positive, and in the 
right half is negative. This means that heat flows from the boundaries towards the 
center of the chain. In Fig. 6(a) we show the heat currents for A = 0.5 and different 
dephasing rates, where the described behavior is observed. The corresponding energy 
profiles are shown in the inset of Fig. 6(a), indicating that the center of the chain is the 
region of lowest energy. In addition, we note that the amplitude of the local currents 
monotonically decreases with 7, and that as 7 increases, so does the amount of lost 
energy, as evidenced from the energy profiles. 

Importantly, the heat current of this configuration satisfies a local diffusion equation 



(J?) = (J xxz > = k X xzV£ xxz , V£f xz = (e i>i+1 ) - (£<_!,<> = (h hl+1 ) - (/^ M >, (26) 



where Ef denotes the local energy density corresponding to XXZ interactions. This 
is verified from our simulations, which indicate that the ratio ( J xxz )/V-E xxz is indeed 
constant in the bulk of the system. Also the obtained conductivity is independent 
of the size of the chain, as expected for a diffusive conductor. For example, for the 
parameters of Fig. 6(a) and 7 = 1, the conductivities obtained for N = 40, 80 and 100 
are Kxxz = 1.997(12), 1.996(7) and 1.996(7), respectively. Therefore the net effect of 
dephasing is to induce diffusive heat transport from the boundaries to the center of the 



Next we study the heat transport through the spin chain when B > and 7 > 0, 



system. 
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Figure 6. Effects of dephasing on the heat transport of the weakly-interacting 
regime, (a) Heat current profile for B = and different dephasing rates. The heat 
conductivities for 7 = 1.00, 0.5 and 0.05 are k X xz = 1.997(12), 3.983(8) and 38(3), 
respectively. The uncertainties correspond to the standard deviation of the values 
of ( J xxz )/Vi? xxz obtained in the bulk of the chain. Inset: Corresponding energy 
profiles, (b) Heat current profile for B = 1 and 7 = 0.5. Both panels correspond to 
N = 40, A = 0.5, / = 0.1 and Y = 1. 



where both (</ xxz ) and (J B ) are non-zero. Recall the characteristic features of each: 
(J B ) = J B is homogeneous through the chain, it has the same properties as the spin 
current and is proportional to the magnetic field, while (</ xxz ) is spatially dependent 
due to dephasing, and is field- independent. In Fig. 6(b) we show the profile of the 
total heat current (J B ) = J B + ( </ xxz ) for a particular set of system parameters. The 
profile is anti-symmetrically centered around J B , and acquires its space dependence 
from (J x ). Near the boundaries, the dephasing processes give energy to the system 
due to the interplay with the driving, so (J B ) < (J^)- Away from the boundaries, 
dephasing induces energy losses from the system, so (J B ) > (J^i). 
A fundamental point to note is that as stated in earlier studies [59, 60], dephasing 
induces a non-equilibrium quantum phase transition from ballistic (7 = 0) to diffusive 
(7 > 0) spin transport in the weakly-interacting regime. Therefore both the spin current 
J s and the heat current J B satisfy a diffusion equation for finite dephasing similar to 
Eq. (25), which reads 

JB = BJ s = KB VEf , VEf = - (&<_!,<) = BVS, (27) 

with Ef denoting the local energy density associated with the magnetic field. Both J s 
and J B are determined by the magnetization imbalance through the system, and have 
the same conductivity, Kb = ks- This leads to a scaling of J B with the size of the system 
like that indicated in Eq. (24). As shown in Fig. 7 this agrees very well with our results. 
The properties of J B and («/ xxz ) indicate that for finite dephasing rates and magnetic 
fields, the total heat current at any site % of the bulk is the sum of two diffusive 
contributions. Each one is determined by the gradient of a different component of 
the local energy density, and has a different conductivity (note the values of Kb and 
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Figure 7. Total output heat current (J B ) ( D ) and its contribution from the magnetic 
field J B (o) as a function of N, for B = 1, / = 0.1, r = 1, A = 0.5 and dephasing 
rates 7 = 0.05,0.50,1.00. The dashed lines are guides to the eye for J-f 1 , while the 
solid lines are fits of the results of J B to the equation J B /AE B = Kb(N — 4)~ Q , 
with AE B = (b N - 2 ,N-i) - (62,3). For 7 = 0.05,0.50,1.00, a = 0.94,1.00,1.00 and 
«b = 15.85, 3.35, 1.86, respectively. Note that the conductivity values differ from those 
of Kxxz, which are field-independent (see Fig. 6). 



Our results indicate that the total heat current does not fulfill a diffusion equation, i.e., 
that an expression such as (J?) = K}iV(Ef+Ef xz ) = hrVE^, with Ef representing the 
total local energy density, is not valid + . Nevertheless, since it consists of the addition 
of two diffusive components, it is clear that the existence of dephasing processes changes 
the nature of the heat transport through the system for A < 1. Therefore, similarly to 
spin transport [59, 60], heat transport shows a dephasing-induced non-equilibrium phase 
transition between ballistic and diffusive behaviors in the weakly-interacting regime. A 
similar result has been reported recently in the non-interacting case A = 0, for a different 
driving scheme and smaller systems [35]. 

Importantly, observe that for any combination of incoherent processes, interaction 
strength and system size, we can always find a magnetic field weak enough to cause 
|(</ xxz )| > J B , and thus that (J$) < 0. Since (Jf) is always positive (both J B > and 
(J_x xz ) > 0), this corresponds to a configuration in which the chain absorbs energy at 
both boundaries and dissipates it by dephasing processes. Some examples of these case 
are shown in Section 5.1. However, in the present work we mostly focus on cases where 
(J B ) > 0, in which the system delivers some energy to the right reservoir absorbed 
from the left reservoir. For such cases we observe the behavior of the total output 

+ This is seen by noting that the ratio ( J/ 1 ) /VEf 1 is not homogeneous through the bulk of the system, 
which results from having kb 7^ «xxz- 



^xxz of Figs. 6 and 7). So we can write 



(Jf ) = J B + ( J xxz ) = n B VEf + k xxz VEJ 



.xxz 



(28) 
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Figure 8. Spin profiles (a) and energy profiles (b) for the same parameters of Fig. 4, 
but with dephasing rate 7 = 0.3. The arrows indicate the increase of the driving /, 
with/ = 0.1, 0.2,..., 1.0. 



current (Jjy); this is shown in Fig. 7 as a function of the size of the system for different 
dephasing rates. Since J B is largely dominant over (J xxz ) for the parameters of Fig. 7, 
the corresponding output heat current («/#) scales with N in a form very similar to J B . 
Finally, it is known that for A < 1 the spin current decreases monotonically with the 
dephasing rate [22, 60, 59]; the same thus happens for J B , as seen in Fig. 7. Since this is 
also the case for ( J xxz ) (see Fig. 6), the total heat current ( J B ) decreases monotonically 
with 7 at any site i of the chain. This shows that dephasing only degrades heat transport 
in the weakly-interacting regime. 

5. Effects of dephasing on heat transport for strong interactions 

In this Section we analyze how dephasing affects heat transport in the strongly- 
interacting regime |A| > 1. We focus on two main effects: the suppression of NDC 
and dephasing-enhanced heat transport. 

5.1. Suppression of NDC 

The NDC effect presented in Fig. 4 is modified by dephasing. This is already made 
clear by considering how dephasing affects the spin and energy profiles, as shown in Fig. 
8 for several drivings and a particular rate 7. The ferromagnetic domains inhibiting 
spin transport at large driving do not appear once 7 > 0, as seen in Fig. 8(a), leading 
to diffusive spin conduction for the entire driving range [22] . The energy profiles, shown 
in Fig. 8(b), feature a more complicated behavior. At driving / < 0.3, their form is 
similar to that of the case without dephasing (see Fig. 4), suggesting heat conduction 
from the left to the right reservoir. At driving / > 0.3, they indicate heat flows from 
both boundaries of the chain towards the bulk. This means that when the driving 
increases beyond a certain value, the output heat current is reversed. 
Now we analyze the heat current. As in the case of weak interactions and 7 > 0, 
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Figure 9. Heat current and its components as a function of the driving /, for 
dephasing rates 7 = 0.05 (o), 7 = 0.30 (□) and 7 = 2.00 (o). The common parameters 
of the simulations are V = 1, B = 0.5, A = 1.5 and N = 16. (a) Total input ((Jf }) 
and output ((Jjy)) heat currents, (b) Components of the output heat current J B and 
(Jjsr )• Note that for a fixed driving, both currents have a non-monotonic behavior 
with 7; for the plotted cases, their largest magnitude occurs for the intermediate 
dephasing rate 7 = 0.3. 

both contributions to the total heat current satisfy a diffusion equation with different 
conductivities, so Eq. (28) remains valid*. We focus on the boundary currents, which 
are shown in Fig. 9(a) for different drivings and dephasing rates. While the input 
current (J^) remains positive and increases its magnitude with /, the output current 
(</Jv) shows a completely different behavior. Initially, for very weak driving, it is positive 
and increases with /. Then it decreases with small increments of the driving. For even 
larger driving, its direction is reversed, so heat starts flowing to the left from the right 
boundary. 

To understand how this behavior of (Jjy) arises, we observe its components J B and 
(J^ xz ) separately, which are shown in Fig. 9(b). The component J B , always flowing to 
the right, features NDC for small dephasing rates (such as 7 = 0.05 for the parameters 
of Fig. 9), but now it presents a finite conductivity at / = 1, a behavior not shown in 
Fig. 9(b). For larger rates (for example, 7 = 0.3 and 2.0), NDC disappears completely 
from J B , its amplitude monotonically increasing with /. In contrast, the component 
(J~n XZ ) always flows to the left, and grows monotonically and rapidly as the driving gets 
stronger. This occurs due to the favored population of high-energy bound states as / 
increases, resulting in a larger amount of energy to be dissipated and thus in a larger heat 
flow from the boundaries to the central part of the chain due to dephasing processes. 
The addition of both components leads to the observed behavior of the total current 
at the right edge of the chain (J B ). Depending on which contribution is dominant, the 
system will feature heat flow from both boundaries towards the bulk, or a net flow to 

* It is known from earlier studies that the spin current shows diffusive behavior for A > 1, at any 
driving with dephasing and at weak driving without dephasing [22, 59]; the same thus applies to J B . 
The diffusivity of {Jf~ ) is proved as before: by observing that the ratio (J xxz )/Vi? xxz is constant 
in the bulk, and independent of N . 
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Figure 10. Dcphasing-enhanced heat transport at weak driving, with / = 0.1, r = 1, 
B = 1 and A = 2.0. (a). Heat current at site N as a function of 7, for several system 
sizes N < 25. (b). Heat currents at the boundaries and their field-dependent and 
independent contributions as a function of 7, for N — 12. Each vertical dashed line 
indicates the optimal dephasing rate 7 op t of the current plotted in the same color; 
7o P t ~ 0.28, 0.39, 0.45, 0.65 for (J$), J B , (jf ), (J$P), respectively. 

the right reservoir. 

5.2. Dephasing- enhanced heat transport 

In earlier work we found that for any driving, the spin current can be significantly 
enhanced by dephasing in the strongly- interacting regime [22]. This result also applies 
to the field-dependent heat current component J B , as seen in Fig. 9(b). Now we 
consider whether a similar effect of environment assistance exists for the total heat 
current in the conventional scenario of unidirectional transport. For different values 
of N, the total output heat current (Jjy) is shown as a function of 7 in Fig. 10(a). 
We find that for chains of small size and moderate dephasing rates, the output current 
can be considerably larger than that of the case 7 = (for example, when N = 8 
the enhancement of the current is ~ 19%), indicating environment-assisted total heat 
transport. For large dephasing rates, the effect is no longer visible due to the degradation 
of coherent dynamics by the frequent perturbation of the environment (e.g. Zeno effect 
[62]) and the bulk energy dissipation. 

The components of the heat current at both boundaries for N = 12 are shown in Fig. 
10(b). We see that in contrast to the weakly-interacting regime (see Fig. 6(a)), the 
boundary currents (Jin) (and in general the currents (Jf ) for all 1) feature initial 
increase and subsequent decay with 7 for strong interactions. So the two components 
of the total heat current show enhancement for moderate dephasing rates. Since both 
flow in opposite directions in the output, they compete with each other to determine 
the total heat delivered to the right reservoir. The corresponding optimal dephasing 
7o P t and the rates of increase and degradation with 7 are different for both; this leads 
to an intricate interplay between J B and (J xxz ) to determine the total output heat 
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Figure 11. Multiple change of monotonicity of the heat transport as a function of 
the dephasing rate 7, with / = 0.1, r = 1, B = 1 and A = 2.0. (a). Heat current at 
site TV as a function of 7, for different system sizes N > 30. (b). Heat currents at the 
boundaries and their field-dependent and independent contributions as a function of 
7, for N = 35. The vertical dashed lines correspond to the dephasing rates 71 s» 0.15 
and 72 ~ 0.6, which indicate the changes on monotonicity for N = 35. 



current. As clearly seen in Fig. 10(b), dephasing-assisted heat transport compared to 
the dephasing-free case exists as long as J B is the dominant contribution to the total 
current over (J^ xz ), and as J B features environment assistance itself. 
Finally, note that in contrast to spin transport, where dephasing assistance compared 
to the case 7 = occurs for any system size [22], the heat transport only displays such 
an enhancement for small chains. For example, for the parameters used in Fig. 10, a 
chain with N = 25 already experiences enough energy dissipation in the bulk to induce 
a barely significant current enhancement (1.3%). Nevertheless, a different assistance 
effect emerges for larger systems, as shown in Fig. 11(a). For small dephasing rates, the 
output heat current ( J B ) is initially degraded compared to the dephasing-free case. But 
from a rate 71, increasing 7 leads to a significant current enhancement, until degradation 
is induced again after a rate 72 due to the Zeno effect. From Fig. 11(b) we can see why 
this multiple change of monotonicity, absent in the input current (J B ), occurs on large 
chains. At small dephasing rates (</ xxz ) is enhanced (negatively) faster than J B , so the 
total current is initially degraded. At 71 the situation is inverted, leading to an increase 
of (J B ) with 7. This tendency continues for moderate dephasing rates up to 72, even 
when both J B and («/ xxz ) are degraded by 7, since the latter decreases faster than the 
former. This shows that a beneficial role of dephasing for heat transport is also present 
in large systems. 



6. Conclusions 



In the present work we have studied the heat transport in non-equilibrium quantum spin 
chains with a homogeneous magnetic field, and observed how it is affected by energy- 
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Figure 12. Schematic diagram of the spin and heat transport regimes; higher color 
intensities correspond to larger currents. In the absence of dephasing (7 = 0), the 
transport is ballistic for the weakly- interacting regime A < 1, diffusive for the strongly- 
interacting regime A > 1 and weak driving, and insulating for A > 1 and strong 
driving. When dephasing is considered (7 > 0), the transport for weak interactions 
becomes diffusive, with the output currents decaying monotonically with 7; this 
corresponds to dephasing degraded transport. In contrast, for strong interactions, 
the transport is enhanced by moderate dephasing rates, until an optimal rate which 
increases with A. As discussed in Ref. [22] for spin transport, the latter behavior 
is expected since as the interactions get stronger, a larger environmental coupling 
is required to break the spin bound states; this results in the boundary shown in the 
figure. A similar qualitative tendency is expected for the heat transport, since the heat 
current and its components are enhanced by moderate dephasing. Further increase of 
7 leads to current degradation, until reaching the Zeno limit, in which the z degrees 
of freedom are frozen and the dynamics of the system is weakly dependent on A. 



dissipating local dephasing processes. The heat current results from a spin imbalance 
at the boundaries of the chain, induced by a Lindblad-type driving, so it corresponds 
to the magnetothermal response of the system. Our results are summarized in Fig. 
12, which qualitatively describes both spin and heat transport. To obtain the steady 
state of the system efficiently, we simulated its time evolution under the corresponding 
Lindblad master equation using the matrix product based TEBD method. We observed 
that the heat current consists of two components: a homogeneous current proportional 
to the magnetic field corresponding to the heat flow carried by the spin current, and 
a field-independent spatially varying current arising from the XXZ kinetics and the 
energy dissipation caused by dephasing, flowing from the boundaries towards the center 
of the chain. 

Initially we found that in the absence of dephasing the heat current is homogeneous 
through the system, and presents the same properties as the spin current. This 
corresponds to ballistic transport in the weakly- interacting regime A < 1, while in 
the strongly-interacting regime A > 1 the system shows diffusive transport for weak 
driving /, NDC as the driving increases, and insulating behavior at large driving / ~ 1. 
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These results hold for any finite amplitude of the magnetic field, thus being independent 
of the (gapless or gapped) nature of the ground state of the Hamiltonian, and instead 
more generally dependent on the energy eigenstructure. 

We proceeded to show that in the presence of dephasing processes, the two components 
of the heat current are finite and compete with each other. Depending on which one 
dominates, the total current at the boundaries can be positive (heat delivered from 
the left to the right reservoir) or of opposite signs (heat flowing towards the bulk). 
We showed that for finite dephasing rates, both current components satisfy a diffusion 
equation with different conductivities and are determined by the gradients of different 
components of the total energy. In the weakly-interacting regime, this corresponds to a 
change of the nature of heat transport, i.e., a dephasing-induced non-equilibrium phase 
transition between ballistic and diffusive regimes. For strong interactions, dephasing 
processes degrade the NDC effect, turning the system into a diffusive conductor for all 
drivings /. 

We also found that for weak interactions, the heat current decreases monotonically as the 
dephasing rate increases, so this coupling to the environment only degrades transport. 
Instead, environment-assisted spin and heat transport emerge at moderate dephasing 
rates in the strongly-interacting regime. While the enhancement of the spin current by 
dephasing exists for any size of the chain [22], that of the heat current only occurs for 
small systems, since dephasing dissipates more energy as the size of the system increases. 
Nevertheless, a beneficial role of dephasing on the total output heat current can also 
be observed for larger systems, when compared to configurations with a small non-zero 
dephasing rate. 

The mechanism underlying the transport enhancement relies on dephasing-induced 
transitions from flat bands of bound states to mobile bands of scattering states [22]. 
Since this is quite generic, we expect that the results reported in this work can be 
experimentally observed for both spin and heat transport, and will be relevant for some 
cuprate materials and controllable interacting systems such as ultracold atomic gases. 
Finally we note that driving schemes alternative to that considered in the present 
work, such as those directly imposing an energy imbalance between the boundaries 
of the chain [3, 55], might result in different transport properties. However, we expect 
dephasing-enhanced transport to emerge also in these cases at strong interactions, given 
the generality of its origin. These effects are currently under study and will be reported 
in a forthcoming publication [68] . 
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